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ABSTRACT 

We develop a semi-analytic model of hierarchical galaxy formation with an improved treatment of 
the evolution of galaxies inside dark matter haloes. We take into account not only dynamical friction 
processes building up the central dominant galaxy, but also binary aggregations of satellite galaxies inside 
a common halo. These deplete small to intermediate mass objects, affecting the slope of the luminosity 
function at its faint end, with significant observable consequences. We model the effect of two-body 
aggregations using the kinetic Smoluchowski equation. This flattens the mass function by an amount 
which depends on the histories of the host haloes as they grow by hierarchical clustering. The description 
of gas cooling, star formation and evolution, and Supernova feedback follows the standard prescriptions 
widely used in semi-analytic modelling. We find that binary aggregations are effective in depleting the 
number of small/intermediate mass galaxies over the redshift range 1 < z < 3, thus flattening the slope 
of the luminosity function at the faint end. At z w the flattening occurs for —20 < Mb < —18, but an 
upturn is obtained at the very faint end for Mb > —16. We compare our predicted luminosity functions 
with those obtained from deep multicolor surveys in the HDF-N, HDF-S, NTT-DF in the rest-frame 
B and UV bands for the redshift ranges < z < 1 and 2.5 < z < 3.5, respectively. The comparison 
shows that the discrepancy of the predictions of other semi-analytic models with the observations is 
considerably reduced at z > 1 and even more at z sa 3 by the effect of binary aggregations. The 
predictions from our dynamical model are discussed and compared with the effects of complementary 
processes (additional starburst recipes, alternative sources of feedback, different mass distribution of the 
dark matter haloes) which may conspire in affecting the shape of the luminosity function. 
Subject headings: galaxies: formation — galaxies: high-redshift — galaxies: interactions — cosmology: 
theory — cosmology: dark matter 

1. INTRODUCTION 



Understanding of galaxy formation has undergone impressive developments in recent years. On the observational side, 
„, ,key step in the study of the statistical properties of high-redshift galaxies has been taken with the discovery of the 
^ Lyman-break" galaxies at z 3 — 4 (Steidel et al. 1996; Adelberger et al. 1998); other important steps have been taken 
with the estimate of the cosmic star formation rate (see Madau, Pozzetti & Dickinson 1998; Fontana et al. 1999), and 
the measurement of the galaxy luminosity functions (LFs) at redshifts up to z sa 4 for UV luminosities extending over 
the range —24 < M <; —18 (Steidel et al. 1999 from spectroscopic data, Pozzetti et al. 1998 and Poli et al. 2001 from 
photometric redshifts) . 

These results, with their unprecedented span in both cosmic time and magnitude, boosted the developement of an 
"ab initio" theory of galaxy formation. The theory has been firmly rooted into the cosmological framework with the 
development of semi-analytic models (SAMs, Kauffmann et al. 1993; Cole et al. 1994; Somerville & Primack 1999; Poli et 
al. 1999; Wu, Fabian & Nulsen 2000; Cole et al. 2000) which link in a single computational structure two main blocks: the 
dynamical history of galaxies as they emerge from primordial dark matter (DM) density perturbations and grow through 
hierarchical merging events; the baryonic processes in the galactic structures, namely, gas cooling, star formation, rise 
and fall of stellar populations, energy feedback from Supernovae. 

Such advances in both theory and observations are pushing the comparison between models and data to higher and 
higher degrees of accuracy. A key role is played by the z-resolved LFs. In fact, the earliest SAMs were calibrated to 
fit the locally observed LF (still with considerable uncertainties in the normalization and in the slope at the faint end), 
and were tested at higher redshift through integrated counts and z-distributions. The recent observations are producing 
a progressive convergence of the local LFs found by different groups (see Zucca et al. 1997, Cross et al. 2001) down 
to Mb ~ —16, with the additional indication that at the faint end the shape is more complex than represented by a 
Schechter form (see Marzke, Huchra & Geller 1994; Loveday 1997). In addition, the resolved LFs now begin to describe 
the evolution of the galaxy population out to z « 4, thus providing a differential test for the model predictions. 

The comparison of the models with such data supports the grand design of hierarchical galaxy formation, but indicates 
that some of the processes included in the SAMs require an improved treatement. For example, the LFs predicted at z « 3 
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overestimate the number of faint (Mi 700 J> —17) galaxies by a factor 5 — 8 when compared with observations based on 
photometric reshifts (see Somerville, Primack & Faber 2001; Poli et al. 2001), and the excess goes beyond the estimated 
incompleteness of the data faintwards of m B <~ 26. A similar trend, though less evident, is found in the redshift range 
around z m 1 (Poli et al. 2001), while at z w the observed upturn in the LF at the faint is not accounted for in the 
current SAM predictions. Addressing the above critical points is crucial to fully understand the physical processes driving 
the galaxy evolution. 

The shape of the predicted LF at faint and intermediate luminosities is affected by two main processes. The first 
concerns the effect of Supernovae (SN hereafter) which heat and partially expell the galactic gas. A stronger SN feedback 
would suppress star formation in smaller haloes thus decreasing their B and UV luminosity and flattening the LF at the 
faint end. However, in the framework of the simple parametrizations currently adopted in SAMs it is extremely difficult to 
increase such feedback without destroying the agreement of the model with other observables; in particular, the predicted 
luminosities for small spiral galaxies would be too faint when compared with existing data concerning the Tully-Fisher 
relation (see Cole et al. 2000). Attempts at non- parametric trcatements of the feedback have been made by Mac Low & 
Ferrara (1999); Goodwin, Pearce & Thomas (2000); Monaco (2001). Different sources of feedback, like that arising from 
the photoionization of the inter-galactic medium by stars and quasars, have been recently included in the SAM framework 
(Benson et al. 2001). 

The other component affecting the shape of the LF is the mass distribution of the galaxies, which is determined by the 
detailed dynamical processes taking place inside the host DM haloes. Among these, tidal stripping and binary aggregations 
of satellite galaxies play a relevant role. However, the former affects mainly the very low-mass end of the mass distribution 
at circular velocities v w 20 — 50 km/s, and has a minor impact on the faint end of the LF (see Benson et al. 2001); 
so we shall focus here on binary aggregations. Such a process is treated in current SAMs only under the approximation 
of orbital decay toward a central dominant galaxy. Here we extend the treatement to include aggregations between all 
galaxies in common DM haloes. On the other hand, we shall adopt the standard SAM prescriptions concerning star 
formation to derive the galaxy LF from the mass distribution; no additional starburst recipes are associated with the 
aggregation events between satellite galaxies. This allows us to single out the dynamical effects of our description without 
introducing new (and uncertain) free parameters in the model, as we discuss in §6. 

The paper is organized as follows: the derivation of the mass function in the framework of SAMs, and the motivation 
for our improved treatment are discussed in §2. A technical description of our approach is presented in §3, where we 
show how we fit our treatment of binary aggregations into the canonical SAM framework. In §4 we briefly recall the basic 
prescriptions that we share with other SAMs to correlate the gas cooling, the star formation and evolution, and the SN 
feedback with the dynamical history of the galaxies. The LFs we derive are compared with the data in §5. In §6 we 
discuss our results and present our conclusions. 

2. THE DYNAMICAL EVOLUTION OF GALAXIES: OVERVIEW AND SPECIFIC MOTIVATIONS 

In SAMs the galaxy mass distribution is derived from the merging histories of the host DM haloes, under the assumption 
that the galaxies contained in each halo coalesce into a central dominant galaxy if their dynamical friction timescale is 
shorter than the halo survival time; the surviving galaxies (commonly referred to as satellite galaxies) retain their identity 
and continue to orbit within the halo. While the histories of the DM condensations rely on a well established framework 
(the extended Press & Schechter theory, EPST, see Bower 1991; Bond et al. 1991; Lacey & Cole 1993), the recipe 
concerning the galaxy fate inside the DM haloes is guided by a posteriori consistency with the outputs of high-resolution 
N-body simulations. 

In reality, additional dynamical processes complement the dynamical friction in driving the evolution of the mass dis- 
tribution. Among these, binary aggregations between satellite galaxies in common haloes have previously been considered 
by Cavaliere, Colafrancesco & Menci (1991, 1992); Cavaliere, & Menci (1993), as a process that would flatten the shape of 
the mass function (and hence of the LF) at small/intermediate masses. This is because such masses aggregate into larger 
units, while the large masses are so few that their binary encounters are unlikely In the above papers, the aggregation- 
driven evolution of the mass distribution was computed in terms of the non-linear Smoluchowski kinetic equation with a 
mass-dependent aggregation rate. The results showed that the effectiveness of the process depends critically on the envi- 
ronment, which in those models was a given input. To describe at the same time the galaxy dynamics and the evolution 
of the host haloes by hierarchical clustering, high-resolution N-body simulations or SAMs are required. 

On the N-body side, recent works (see Klypin et al. 1999 for pure DM, and Murali et al. 2001 for hydrodynamical 
simulations) indicate a complex galaxy growth. While at very low circular velocities (v w 20 — 50 km/s) tidal stripping 
may affect the mass distribution (Gnedin & Ostriker 1997; see also Taylor & Babul 2000; Taffoni et al. 2001), at larger 
masses (v ~ 100 km/s) binary aggregations of satellite galaxies do play a relevant role which complements the coalescence 
into a central galaxy through dynamical friction. Indeed, Murali et al. (2001) analyze the competing role of the two 
growth modes of the simulated galaxies in terms of an evolutionary equation which includes also a binary aggregation 
term of the Smoluchowski type. 

On the SAM side, efforts to insert the aggregations between satellite galaxies have recently been started by Somerville 
& Primack (1999). However, they adopt a cross section (derived from the N-body simulations of galaxy encounters by 
Makino & Hut 1997) valid for equal galaxies in clusters with velocity dispersion much higher than the internal galaxy 
dispersion, a condition that allowed them to adopt a one- body treatement for the aggregations. 

Here we take a step forward, and consider in closer detail the actual two-body dynamics of aggregations; we shall also 
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adopt a cross section valid down to small groups with velocity dispersions close to those internal to galaxies, as is the case 
at early times in the hierarchical clustering picture. 

In fact, we develop a SAM including both dynamical friction and binary aggregations. Instead of employing a Monte 
Carlo simulation, as usual for SAMs, we follow the evolution of the galaxy mass distribution by solving numerically a 
set of evolutionary equations, as in Poli ct al. (1999). The subset describing the two-body dynamics is constituted by 
the (non-linear) kinetic Smoluchowski equation. This modifies the mass function resulting from dynamical friction by an 
amount which depends on the properties of the host DM haloes, which in turn evolve according to the EPST theory. The 
basic description of gas cooling, star formation and evolution, and SN feedback is kept unchanged with respect to the 
standard SAM prescriptions given, e.g., by Poli et al. (1999) and Cole et al. (2000). 

We first compare our results with those from existing N-body simulations, and then compare our LFs with those 
obtained from spectroscopic and deep imaging surveys down to faint magnitutes {Iab ~ 27.2), at redshifts up to z ~ 4. 
This allows us to test the effects of binary aggregations on the SAM predictions over a wide range of cosmic times and 
masses. 

3. THE GROWTH OF THE GALAXY MASSES 

We consider the number N(m, M, t) dm dM per Mpc 3 of galaxies with mass in the range dm about m, residing in haloes 
with masses in the range M to M+dM at cosmic time t. At the initial time to we assign one galaxy to each halo, a condition 
which formally translates into N(m, M, to) = Nh(M, to)8(m — M) where 8 is the Dirac delta function. Our default choice 
for the halo mass distribution Nu(M,t) is the standard Press & Schechter (1974) expression, whose dependence on the 
cosmological parameters and on the spectrum of primordial density perturbation is recalled in Appendix A; but we shall 
also explore the effects of adopting the Sheth & Tormen (1999) mass distribution, also recalled in Appendix A. 

Given the merging history of DM haloes described by EPST, two main processes affect the evolution of N(m, M, t): the 
orbital decay of satellite galaxies onto a central dominant galaxy due to dynamical friction, and the binary aggregations 
between satellite galaxies. In the following two subsections we investigate the effects of the above processes on the evolution 
of N(m, M, t). A schematic view of the two dynamical processes driving the growth of galaxies in our model is given in 
fig. 1. The corresponding evolution of the mass distribution of galaxies is presented in the last subsection 3.4. 




Fig. 1. - A schematic representation of the various terms contributing to the evolution of the galaxy velocity function N(v,V) inside DM 
haloes of circular velocity V . The first three processes marked with DF correspond to the construction and destruction terms (eq. 1) from 
dynamical friction following the merging of the host haloes (see text). The last two processes, marked with BA, represent binary aggregations 
of satellite galaxies (see eq. 3) inside the DM host halo. Their combined effects drive the evolution of N(v, V,t). 
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3.1. Coalescence driven by dynamical friction 

In the following we shall frequently adopt the circular velocities as the proper quantity to mark the depth of the potential 
wells of galaxies and of their host haloes. As for the latter, this is related to the mass through the relation V = \J GM/R 
where the limiting radius R is the radius within which the mean mass density is 200 p c , where p c is the critical density 
at the redhsift z where the halo is identified. The relation between R and M is given, e.g., in Navarro, Frcnk & White 
(1997) and Mo, Mao & White (1998) and takes the form R = 1.63 10~ 2 (M/h- 1 Mo) 1 / 3 [n /^(z)]~ 1/3 (1 + z)- 1 h- 1 kpc 
for Einstein-de-Sitter (fl = 1, = 0), open (fl < 1, = 0), and flat (CIq + fl\ = 1) universes. 

As to a galaxy inside a host DM halo, the circular velocity v is related to m through m — v 2 rud/G; here rud is the tidal 
radius, within which the mean density of the galaxy exceeds the average density of the host halo interior the pericentre of 
the galactic orbit. This ensures that the galactic subhalo survives the tidal stripping due to the host halo potential wells 
and retains its own identity. Our computation of rud is given in Appendix B. 

The above relations allow us to relate the mass and the circular velocity distributions by applying the proper Jacobian. 

The evolution of N(v, V, t) driven by dynamical friction is computed as follows. At each time step, we first compute 
the conditional probability d 2 Pff(V' — ► V, t)/dV dt that a given halo with circular velocity V at time t has a progenitor 
with circular velocity V at t' = t — At; this is calculated in the framework of the EPST, and its expression (depending 
on cosmology and on the perturbation spectrum as given in Bower 1991; Bond et al. 1991; White & Frcnk 1991, Lacey 
& Cole 1993), is recalled in Appendix A for a generic t' (see eq. A-3). Similarly, from the EPST we compute the inverse 
conditional probability d 2 Pn{y — > V ,t)/dV dt that a given halo with circular velocity V at time t ends up in a halo 
with circular velocity V at t + At (see eq. A-2). 

According to the canonical prescriptions usually adopted by SAMs, during halo merging a galaxy contained in one of 
the parent halos contributes to enrich the dominant galaxy of the common halo of circular velocity V if its coalescence 
time Tdf(v) is shorter than the halo survival time n{V) predicted in the EPST. The circular velocity of the merger is then 
equated to that of the host halo. Or else, for Tdf > n, the galaxy retains its identity. Thus the increment in the number 
N(v, V, t) of galaxies of given v in halos of given V is linked to the increment of DM halos with V, whose contribution from 
progenitors V' is Nh(V, t) d 2 Ph{V — > V, t)/dV'dt; the link occurs through the probabilities: i) of forming one dominant 
galaxy of circular velocity v = V by coalescing galaxies contained in the parent halos with lower velocities v' < v; ii) of 
finding galaxies with velocities v which have not coalesced despite the merging of their host haloes. The corresponding 
decrement is due to the inclusion of galaxies with current velocity v into a larger halo V > V. The construction and 
destruction terms described above are schematically represented in the upper part of fig. 1 . The corresponding evolution 
of N(v, V, t) in a single timestep At is expressed by 

N(v,V,t + At)-N(v,V,t) = 
= AtS(v V) jT dv' £ dV'N H (V',t) dPHi Jy^ t V,t) ^^- P rob[r df (v') < n(V)] 

+ * [ dv > itv>, {1 _ M) < Tl{v)] } 

where Nt(V) — J dv' N (v' , V, t) is the total number of galaxies per Mpc 3 , per unit halo rotational velocity. The first 
term on the right-hand side implies that if any galaxy with circular velocity v' < V is included in a halo with circular 
velocity V and its orbit decays (due to dynamical friction) to the center of such a halo within the halo survival time, then 
a central galaxy with velocity v = V is formed inside the halo. 

Coalescence is caused by loss of galaxy energy and orbital angular momentum due to dynamical friction to the halo 
material. The timescale of such process is usually determined through orbit averaging of the Chandrasekhar formula (see 
Lacey & Cole 1993); the result is given in Cole et al. (2000) as follows: 

0.3722 M 

Tdf = O T dyn . (2) 

ln(M/m) m 

Here Td yn = n R/V is the dynamical time of the halo, and O = [J/ J C (E)} - 78 [r c (R)/R] 2 contains the dependence on the 
initial energy E and angular momentum J of the galaxies, in terms of the angular momentum J c and the radius r c of 
the circular orbits corresponding to E. Although the values of are statistically distributed (approximately following a 
lognormal function with (logio O) = —0.14 and ((logio O — {logw O)) 2 ) 1 / 2 w 0.26, see Tormen 1997), in the following we 
shall use for O its average value. Note that the orbit averaging procedure (adopted to derive eq. 2) is not appropriate for 
very elongated orbits for which the local density varies on a much shorter timescale than that given by eq. (2); however, 
such orbits are expected to be a minority (<J 15%, see Ghigna et al. 1998). 

The probability prob[r < tj(V)1 for a DM halo of velocity V to have a survival time r; larger than a given value r has 
been computed by Lacey & Cole (1993, see their eq. 2.21) in the framework of the EPST, and is recalled in Appendix A. 
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3.2. Binary aggregations of satellite galaxies 

In addition to the above coalescence process, we include binary aggregations. So for each DM halo with circular velocity 
V, we compute the further evolution of the galaxy mass distribution in a timestep due to binary aggregations. This is 
described by the Smoluchowski equation, which we write in terms of masses for the sake of simplicity: 

N(m, M, t + At)- N{m, M, t) = 
1 t m 

-At dm' N(Tn',M,t)N(m-m',M,t)T~ 1 (m' 7 m-m',V) 



- AtN{m,M,t) / dm' N(m\M,t)T-J- (m,m',V) , (3) 



T agg 





where T~ gg {m,m' ,V) is the aggregation rate depending on the DM halo where the galaxies m and m! reside. The first 
term describes the construction of galaxies with mass m from smaller ones with mass m! and m — m', while the second 
represents the destruction of galaxies m due to their aggregation with others. A schematic representation of the two terms 
governing the binary aggregations is given in the lower part of fig. 1. 

The aggregation rate is governed by r~ gg = £ V re i/ (AitR 3 /?>) (see Cavaliere, Colafrancesco & Menci 1992). Here V re i is 
the average relative velocity of galaxies in the DM halo whose rms value is equal to twice the halo 1-D velocity dispersion 
ay ~ V I y/2, and £ is the cross section. The latter is given for nearly grazing, weakly hyperbolic encounters by Saslaw 
(1985), and by Cavaliere, Colafrancesco, & Menci (1992). It includes a geometrical term (proportional to the area of the 
galaxies), and a focussing factor ~ 1 + {v/V re i) 2 that accounts for the enhancement of £ in slow encounters with resonance 
between the internal and the orbital degrees of freedom (see also Binney & Tremaine 1986). Thus, the average rate for 
binary aggregations is 

where the average is over the relative velocities V re i. The distribution of the encounter velocities is assumed to be 
Maxwellian, namely 

g (y rel ) = ./? Yk e- v ^v . (5) 

Note that for encounters between equal galaxies with r' — r and N(M') — N(M) = N the scaling of the aggregation rate 
(||) reduces to T~ gg oc (r 2 (1 + v 2 /V 2 el ) V re i). Performing the average over the distribution g(V re i) in terms of the rescaled 
variable y = V re i/v yields r~ gg oc r 2 v 4 ay 3 R~ 3 1(x), where the function I(x) = J dyy 3 exp(—y 2 /x 2 ) + f dyy exp(—y 2 /x 2 ) 
tends to a constant when the ratio x = cry/v — > oo. Thus the aggregation rate (^) reduces to the expression given by 
Makino & Hut (1997) - originally derived analytically by Mamon (1992) - and the r.h.s. of eq. (3) becomes proportional 
to iV 2 r 2 v 4 ay 3 R~ 3 (the effective rate adopted by Somerville & Primack 1999) in the proper limits of large encounter 
velocities relative to the internal galaxy velocity dispersion, and of encounters between equal galaxies. 

Finally, note that after performing the average of eq. (4) over the distribution g(V re i) only ay enters the computation; in 
other words, we use an average description for encounter velocities typical of the environment considered. The geometrical 
cross section contained in eq. (4) does not properly describe single events with V re i S> v, but these are rare for our typical 
values of ay /v in the range from 1 to about 4. In cases with larger ratios ay jv the effect of the aggregations is suppressed 
since T agg oc a v as shown above, but even inside rich clusters the averaged cross section that we adopt is consistent with 
the simulation results, see Makino & Hut (1997); for numerical computations of the role of interactions inside clusters see 
also Lanzoni (2000). 

In sum, our cross section provides an accurate average description of aggregations in systems where the velocity disper- 
sion is close to the galaxy circular velocity and the focussing term is relevant; on the other hand, it constitutes an effective 
average approximation for encounters in environments with large velocity dispersion up to the scale of rich clusters, where 
the aggregations are disfavoured anyway. As a final global check, we have verified that the insertion of an artificial cutoff 
at V re i — 4 v in E does not change our results. 

3.3. Numerical solutions of the equations: test cases 

The equations (0) and (||) describing the evolution of N(m, M, t) are integrated numerically on a grid of circular 
velocities and cosmic times with step At = 10~ 2 H^ 1 . 

In order to test our numerical code, we first run the computation for two relevant simple cases where analytic solutions 
are available. In the limit tm — > (i.e., when merging of the host haloes is promptly followed by coalescence of the 
galaxies within them by dynamical friction) the solution N(v, V, t) of eq. (1) when integrated over the circular velocity V 
must yield the Press & Schechter mass distribution. The comparison between the numerical and the analytic solution in 
this case is performed at three different times in the top panel of fig. 2, which shows that the numerical solutions remain 
close to the analytic Press & Schechter form over the whole range of v; in fact, the relative deviation is always smaller 
than 5%. 

To test the code section concerning the Smoluchowski equation, we numerically solve eq. (3) in the case of constant 
aggregation rate T~ gg = const for galaxies within a host halo of given mass M. In this case the exact solution (Smoluchowski 
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1916; Trubnikov 1971) has the form N(m,t) — [A /m~.(t)] exp[—m/m*(t)], where the constant A a is related to the total 
mass M. contained in the system, and m*(i) = m*o if /to) is a characteristic mass linearly growing with time. Trubnikov 
(1971) has shown that such a solution holds for very general initial conditions after a transient time. The comparison 
with the numerical solution is performed in the bottom panel of fig. 2, where the initial condition (dotted line) has been 
chosen to be N(m,to) oc (m/m^o)" 1 exp[—m/m*o] normalized as to yield a total mass M. = 1.54 10 2 m»o. Note how the 
numerical solution progressively flattens at small masses during the transient, to approach the exact solution (computed 
at t = 3 to) with a relative error smaller than 3% over the whole mass range. Note also that the numerical solution 
conserves the total galaxy mass with high accuracy. 
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Fig. 2 - Comparison between the numerical and the analytic solutions of eqs. (1) and (3) for two test cases. Top panel: the numerical 
solutions of eq. (1) in the limit of vanishing Tjf (dashed curves) are compared with the Press &; Schechter form (solid curves). The curves 
refer to redshifts 2 = 8, 2.6 and (from left to right at the high-mass end). Bottom panel: the numerical solution of the Smoluchowski eq. 
(3) , computed for a constant aggregation time (dot- dashed curves) at 6 equally spaced time intervals from to to 3 to , are compared with the 
corresponding analytic solution (given in the text) computed at t = 3 to (heavy solid line). The dotted line represents the initial condition. 
The inset shows the variation with time of the total mass M corresponding to the numerical solution. 

3.4. The evolution of the galaxy circular velocity distribution 

Having tested our code, we proceed to compute the complete evolution of N{v, V, t). At each time step, we first compute 
the change of N(v, V, t) due to dynamical friction (the right hand side of eq. |l|), and then the further change due to binary 
aggregations of satellites (the right hand side of eq. ||) using the physical values of r^/ and T agg given in eqs. (2) and (4). 
A discussion of the role of the two terms is given in the final §6. 

Once the complete evolution of N(v,V,i) is found for all times, we compute the probability of finding a galaxy with 
circular velocity v in a halo of circular velocity V as f(v,V,t)dv = N(v,V,t)dv/NH(V,t). Then the total density 
distribution of galaxies with circular velocity v is computed from N(v,t) — J dV Njj(V,i) f(v,V,t), where the halo 
distribution Nh(V, t) takes the canonical Press & Schechter form. This is the distribution of circular velocities irrespective 
of the halo to which the galaxies belong (the global velocity distribution) . 

The "transition probability" for galaxies is given by p(v', t', v, t) = f™ dV dV ^XL^YS. V ' , t') f(v, V, t), 
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where ifsiY-^^jXiH [ s the fraction of mass in haloes of circular velocity V at time t', and later at t > t' in haloes with 
V; this is given by EPST (see Appendix A). 

The evolution of the galaxy circular velocity distribution N(y,t) resulting from the full dynamics is shown in fig. 3 
for a CDM power spectrum of primordial density perturbations and for our reference set of cosmological/cosmogonical 
parameters: £Iq = 0.3, fl\ = 0.7, and Hubble constant h = 0.7 in units of 100 km/s/Mpc. For comparison, we also show 
the evolution corresponding to dynamical friction alone (as usually considered in SAMs), and the evolution of the halo 
velocity distribution as given by the Press & Schechter formula. 




1 10 1 10 

v(100km/s) v(100km/s) 



Fig. 3 - The galaxy velocity function at four different redshifts. The velocity is measured in units of 100 km/s. The solid line is the velocity 
function resulting from the full model including dynamical friction and satellite aggregation (eqs. 1 and 3). The function N(v) resulting 
from dynamical friction alone is shown as a dashed line, while the velocity function of the host haloes (computed from the Press & Schechter 
formula) is also shown for comparison as a dotted line. In the bottom-right panel we also show the velocity function resulting from the N-body 
simulations by Klypin et al. (1999); the solid squares correspond to a simulation with a box size of 60 h~ 1 Mpc, and the open circles to a box 
size of 7.5 h~ 1 Mpc which allows for a finer mass resolution. Both cases are computed for A-CDM cosmology with parameters as in the text. 
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First consider the effect of dynamical friction. Compared to the halo velocity distributions, the galaxy distribution 
shows a delayed evolution for z < 3, particulary evident at low rcdshift. This is due to the longer timescale on which 
dynamical friction occurs (eq. ||) as compared to the survival time of the halo (which is of order Td yn )- As a consequence, 
while at high rcdshift the merging of haloes is promptly followed by coalescence of the galaxies within them, at later 
times the number of galaxies accumulating in the haloes increases due to the longer time taken by dynamical friction to 
occur. At small masses this implies more galaxies than haloes, while at large v the increase in the number of haloes is 
not followed by a corresponding increase of massive galaxies since these are formed on a longer timescale. 

The effect of binary aggregations is to flatten the slope of the galaxy velocity distribution at the low-mass end (v ~ 
50 — 150 km/s). Such a process is already effective at z ~ 2.5 and continues down to z « 0, where our N(v) is fully 
consistent with available results from N-body simulations as shown by fig. 3. However, in the range from z ~ 1 to z « 0, 
the aggregations deplete more efficiently objects with intermediate mass (v ~ 100 — 150 km/s), thus producing an upturn 
of the local velocity distributions at low v. This is because at low z the galaxies are hosted in haloes that are typically 
more massive than at higher z, according to the hierarchical clustering. The large relative velocities V re i that galaxies 
have in such deeper potential wells suppress the aggregation of galaxies with low v due to the v/V re i term in the cross 
sections (||), as stressed by Cavaliere & Menci (1997). 

To further assess the consistency of our results with N-body simulations and to guide our interpretation of the aggre- 
gation role we plot in fig. 4 the overall merging rates predicted by our model and compare them with those obtained 
in recent hydrodynamical simulations in a cosmological framework (Murali et al. 2001). In the top panel we show the 
net number increment as a function of the cosmic time for all galaxies with baryonic mass larger than 5.610 10 Mq (the 
baryonic mass associated with the DM mass in our model will be derived in §4). We also show as a dashed line the 
construction rate (i.e., the sum of the positive terms in eq. |l| and jj), also illustrated in fig. 1), to show that it dominates 
over the destruction terms for such galactic masses. 
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Fig. 4 - Contributions to the number density of galaxies with masses above 5.6 10 10 Mq. Top panel: the net rate (creation minus destruction) 
predicted by the full model including dynamical friction and aggregations of satellite galaxies is shown as a solid line. This is compared to 
that obtained from the simulations by Murali et al. (2001, represented by the dots) with the same cosmological parameters adopted here. 
The dashed line shows the construction rate alone, as predicted by the model. Bottom panel: the destruction rate (the difference between the 
two curves in the top panel) resulting from the model (solid heavy line) is compared with the simulations (dots). The contributions to such 
destruction rate from dynamical friction (dotted line) and from aggregations of satellite (dashed line) are also shown. 
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Notice how simulations and our model agree in yielding a net rate which declines with time after a peak at z as 2.5 — 3. 
The decline of such a rate for z < 1 is consistent with the observed one (see LeFevre et al. 2000; Carlberg et al. 2000) 
within the uncertainties entering the comparison. 

To assess the relative role of dynamical friction and satellite aggregations we show in the bottom panel the destruction 
rate for the same mass threshold, compared with the N-body results. Note the smooth decline at low redshift in both 
the simulations and the model (the heavy solid line) after the peak at z as 2.5. Fig. 4 shows in detail that the position 
of the peak is mainly determined by dynamical friction; satellite aggregations begin to contribute significantly to the 
destruction rate at z as 2.5 — 3. Then they sustain the destruction rate during its slow decline to lower z. Thus, first the 
hierarchical clustering drives the construction (top panel) and the destruction (bottom panel) of the basic galactic blocks 
within a short stretch of cosmic time. Then binary aggregations further deplete the number of galaxies. The evolution 
of the destruction rate corresponding to binary aggregations is due to the shifting balance of different terms: first, the 
aggregation efficiency increases due to the growth of the galaxy sizes and of the number of galaxies contained in the host 
halos (see the aggregation rate in eq. ^); at later times, the above terms are balanced and/or overwhelmed by the increase 
of the galaxy relative velocities, which grow in time following the increasing mass of the host haloes and tend to suppress 
the binary aggregation rate eq. (|j). 

Such a picture is confirmed by the fact that when a lower mass threshold is considered, the high-z shape of the 
destruction rate remains similar to that in fig. 4, while the low-z tail is appreciably suppressed. This is because for small 
galaxies the effect of the increase of the relative velocities with time is stronger, and prevents binary aggregations from 
sustaining the destruction rate at small z. So at z as the velocity functions in fig. 3 has an upturn at small masses. 

Having discussed the full dynamics, we now turn to recall the basic points of the stellar section of our SAM. 

4. STAR FORMATION AND EVOLUTION 

To link the stellar section of the model with the dynamics we adopt the standard procedure commonly used in SAMs 
(Somerville & Primack 1999; Cole et al. 2000); we shall give here a brief presentation of how such a procedure fits into 
our statistical treatement of the evolution of the mass distribution of galaxies inside host haloes. 

The baryonic content (J7fc/r2 m )m of the galaxy is divided into a hot phase with mass m/, at the virial temperature 
T = (1/2) limn v 2 Ik {ma is the proton mass and /i is the mean molecular weight), a cold phase with mass m c able to 
radiatively cool within the galaxy survival time, and the stars (with total mass m*) forming from the cold phase on a 
time scale r*. 

Initially, all the baryons (for which we adopt a density parameter Vt^ = 0.02) are assigned to the hot phase. Then, for 
each galaxy mass m (with circular velocity v and radius r) we compute the baryons in the three phases as follows. 

a) The cold phase. Its mass m c is increased from inside out by cooling processes: at each time step At, we have then 
Am c {v) — 4ir r 2 ool p g Ar coo i; for the gas density we adopt the form p g cx l/(r 2 + r 2 ore ), with the value of r core determined 
by requiring that the density at the virial radius is the same that would have been obtained has no gas cooled (hence 
depending on the past history of m c and corresponding to each galaxy mass m at time t). The cooling radius r coo ; 
at time t is computed by equating the cooling time r coo / = (3/2) kT '/ fimn pg(f) A(T) to the current time t. The cooling 
function A(T) is taken from Sutherland & Dopita (1993) for a mixture of 77 % hydrogen and 23 % helium. 

The gas which cools settles into a rotationally supported disk; following Mo, Mao & White (1998), its radius rd(v) 
and the rotation velocity Vd(v) are related to the galaxy circular velocity assuming that the angular momentum of the 
baryons which settle into the disc is a fixed fractions jd of the total angular momentum J of the galaxy dark halo. The 
latter is usually expressed in the adimensional form A = J j(\E\^l 2 jGrn^l 2 ). We use the relations Td(v, A, c, m<2, jd) and 
Vd(v, A, c,md, jd) given by the above authors for a fraction of cold gas = m c /m in a Navarro, Frenk & White (1997) 
potential with concentration c(v). In the following, such relations will be used only to compute the Tully-Fisher relation 
resulting from our model (see fig. 5 below) and the dust extinction to be applied to the galaxy integrated stellar emission 
(see eq. [7] and below). In such computations we will assume jd — 0.05 (as discussed by the above authors) and integrate 
over the full distribution of A (a lognorm expression with mean A = 0.05, and dispersion <j\ = 0.5 in log(X), see Warren 
et al. 1992; Cole &: Lacey 1996; Steinmetz & Bartelmann 1995) to obtain average values for Vd and r^. The disk sizes 
resulting from the above values of jd and are consistent with recent observations, as shown in Giallongo et al. (2000). 

b) The stars: the mass of baryons turned into stars in the timestep At is assumed to be Am* = Aim c /r*, with a 
star formation timescale r* = e^ 1 (iv/200km/s) a * proportional to the dynamical time of the disk Td = rd/vd- The 
normalization e* and the exponent a* are free parameters. Our default choice for them is e* = 0.05 and a* = 1.5, as in 
the fiducial model of Cole et al. (2000). This allows us to directly compare ours with previous results, and to single out 
the effects of the dynamics on the LFs. Note that the star formation is linked only to the changes of the cold material 
following the dynamical history of the galaxies; as discussed in §6 we do not include possible starbursts associated with 
newly formed mergers. 

c) The hot phase is replenished by part of the cool baryons, namely, those re-heated and ejected into the hot phase by 
SN feedback. This is the most uncertain among the baryonic processes included in SAMs, and is currently parametrized 
by the simple expression Amn = Am* (v/vh) ah relating the amount Amu of the mass re- heated in a time step At to the 
mass of stars Am* formed. The exponent ah, a free parameter, models the increased feedback efficiency in galaxies with 
decreasing mass. The choice of at has a considerable impact on the faint end of the LF (larger values yielding flatter 
shapes), but its values are strongly constrained by the observed Tully-Fisher relation. In fact, increasing ah produces 
fainter luminosities for small mass objects. This flattens the resulting LF at the faint end, but it also moves the predicted 
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Tully-Fisher relation at small v away from the observed range. Our default choice of the parameters is ah = 2 and 
Vh = 200 km s -1 , since these yield the best joint fit to both the local LF and the Tully- Fisher relation, as shown below. 

After having computed Ato c , Am, and Am/, from the processes a), b), c), we compute the increments due to galaxy 
coalescence/aggregation. The probability p(v' , t', v, t) that a galaxy with circular velocity v' at time t' is included into a 
galaxy with velocity v at time t has been calculated in §3.4. Then, the average star content in a halo with circular velocity 
v is updated over the time grid according to following equation: 

m*(t>,i) = y"/ dv' N J\ ' ^ p(v', U,v,t) Am«(i)', U) , (6) 
• Jo 1 ^\ v , t ) 

where the sum extends over the time steps i$ = to + i ■ At ranging from the initial time to to t. Analogous equations 
govern the increments of m c and mt- Such a procedure is similar to that introduced by Frenk & White (1991); however, 
here the transition probability for the galaxies is computed on the basis of the dynamics described in §2. 

From the latter equation, the corresponding integrated stellar emission S\(v,t) at the wavelength A is computed by 
convolving with the spectral energy distribution (f>\ obtained from population synthesis models: 

S x {v,t)= f dt'<j> x (t-t')m4v,t') . (7) 
Jo 

Note that the average star formation m*(v,t') of galaxies with circular velocity v at t' is that corresponding to the star 
mass computed in eq. (6). Substituting its expression into eq. (7) demonstrates that S\(v,t) contains the integrated 
contributions of star formation in the smaller progenitor systems (with circular velocity v' at times t' < t) which by 
the time t have been included in the haloes with circular velocity v. The integrations over the time t' and the velocity 
v' entering eq. (7) account for the average building up of the stellar population in hierarchically growing galaxies, as 
described by Frenk & White (1991). 

In the following we adopt <p\ taken from Bruzual & Chariot (1993), with a Salpctcr IMF. The dust extinction affecting 
the above luminosities is computed assuming the dust optical depth to be proportional to the metallicity Z co id of the 
cold phase (computed assuming a constant effective yield and that the metals are re-ejected to the hot phase in the same 
proportion as the reheated gas Ami) and to the disk surface density, so that for the V- band Ty oc m c Z co id/irrj. The 
proportionality constant is taken as a free parameter chosen as to fit the bright end of the local LF (see below and fig. 
5); this yields for the proportionality constant the value 3.5Mg 1 pc 2 when the stellar yield is such as to produce a solar 
metallicity for a v = 220 km/s galaxy. To compute the extinction in the other bands different extinction curves will be 
considered, including the Milky Way (MW), the Small Magellanic Cloud (SMC), and the Calzetti extinctions, see Calzetti 
(1997). 

-10 I 1 I 1 1 1 I 1 1 1 I 1 1 1 I 1 1 1 I ' I -10 I 1 I 1 1 1 I 1 1 1 I 1 1 1 I 1 1 1 I ' 
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Fig. 5 - The Tully-Fisher relation relating the disk rotation velocity to the luminosity of spiral galaxies. Left panel refers to the our fiducial 
case = 2, while the right panel to the high-feedback case = 5. The shaded areas represent the region of the Mi — v plane allowed by the 
observations (Mathewson et al. 1992; Willick et al. 1996; Giovanclli et al. 1997). The solid curves represent the model predictions assuming 
the disk rotation velocity equal to the DM circular velocity v, while the dashed curves are computed adopting the disk circular velocity v^(v) 
as a measure of the disk rotation velocity. 
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The above equation (7) relates the observable properties of galaxies (specifically their integrated stellar emission S\) to 
their dynamical properties. The luminosity-circular velocity relation constitutes a first, key prediction that can be tested 
against the observed Tully-Fisher relation. The comparison between the outcomes of our model (with our default choices 
for the free parameters) and the observations is shown in the top panel of fig. 5. The agreement is satisfactory, although 
the predicted magnitudes at given v are slighly fainter than the data when the disk circular velocity Vd(v) is used as a 
measure of the disk rotation velocity, as appropriate. 

To show that increasing ah is not a viable way to get flatter LFs (as anticipated above at point c), we also show the 
Tully-Fisher relation with ah = 5, the value initially adopted by Cole et al. (1994) to get a flat local LF; in this case the 
predictions at faint luminosities fail to match the observations. Other 

5. THE EVOLUTION OF THE GALAXY LUMINOSITY FUNCTION 

From the galaxy velocity distribution N(v,t) derived in §3 (see fig. 3) and the S\ — v relation discussed above we 
derive the galaxy LF predicted by our model. The results are compared with the observational LF obtained from deep 
surveys. In particular, our predictions at high z are compared with the observational LFs derived by Poli et al. (2001) 
using photometric redshifts for galaxies in the ESO New Technology Telescope deep field and in the Hubble Deep Field 
North and South, down to the limiting magnitudes Iab = 27.2. The LFs have been computed in the rest-frame B band 
for < z < 1, and at the rest- frame wavelength of 1700 A for higher z. We stress that in the magnitude range where 
the photometric data overlap with those from spectroscopic surveys (see, e.g., Steidel et al. 1999), the photometric and 
spectroscopic LFs agree to a remarkable degree of precision. 

5.1. Comparison with data: the effect of satellite aggregations 

The comparison between the observed and the predicted LF is given in figs. 6 and 7. We compare with the data both 
the LFs correponding to the standard model (where only coalescence driven by dynamical friction is present) and those 
from the complete dynamics including binary aggregations between satellite galaxies; hereafter we shall refer to such cases 
as DF and DF+BA. In both cases, we adopt our reference set of cosmological parameters (given in §3.4) for a Universe 
dominated by a cosmological constant. 
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Fig. 6 - The local LF computed with the DF process alone (left panel) and with the inclusion of satellite aggregations (DF+BA, right panel). 
The different curves (almost overlapping here) refer to different choices of the dust extinction curve; MW (dashed line), SMC (dot-long dash), 
and Calzetti (dotted). The data are taken from Zucca et al. (1997, filled squares) and from the 2dFGRS survey (Madgwick et al. 2001, open 
circles). 



The low-z LFs presented in fig. 6 show how satellite aggregations flatten the LF at faint/intermediate luminosities. 
Such an effect is purely dynamical, so that the agreement of the model with the Tully-Fisher relation shown in fig. 5 is 
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unchanged. At the very faint end, note the upturn of the predicted local LF in the DF+BA case, as expected from our 
discussion in §3.4. The local trend toward a flatter LF shown by the DF+BA model persists at z ~ 1 (see fig. 7), where 
our model actually provides a better fit to the shape of the data distribution when compared to the DF case. 

The high-z LFs are compared with data in the bottom panels of fig. 7. Note that all models overpredict the number of 
faint galaxies as noticed in Poli et al. (2001). Several processes not yet properly inserted in the SAMs can be at the origin 
of the discrepancy. As discussed in sect. 4, simple variations of the stellar feedback prescritions adopted in the model 
are not a viable solution, being constrained by the observed Tully-Fisher relation. On the other hand, our treatement of 
binary aggregations considerably reduces the discrepancy, as shown by the right-hand panels in fig. 7. A residual excess 
of the predicted over the observed LF remains at the faint end. 




Fig. 7 - The LF in the B-band at z = 0.7 (top panels) and at 1700 A at z = 3 (bottom panels) for the DF process alone (left panels) and with 
the inclusion of satellite aggregations (DF+BA, right panels). The solid curves refer to the LF without dust extinction, while the other curves 
refer to different choices of the dust extinction curve: MW (dashed line), SMC (dot-long dash), and Calzetti (dotted). The data with the error 
bars correspond to the LF computed from photometric redshifts by Poli et al. (2001); the data corresponding to the spectroscopic survey by 
Steidel et al. (1999) are also shown by the circles in the bottom panels. In the right-bottom panel we also show the effect of adopting the Sheth 
& Tormen (1999) instead of the Press & Schechter expression for the halo mass distribution; the corresponding galaxy LFs (derived including 
DF+BA) are represented as thin curves (the line types correspond to the different extinction curves as above). 
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A possible origin for it can be found in the statistics of the DM haloes (the Press & Schechter mass distribution) hosting 
the galaxies. To investigate the issue we computed the LF from our model (including the binary aggregations of satellite 
galaxies) adopting the Shcth & Tormen (1999) form (recalled in the Appendix A) for the mass function of the DM haloes 
hosting the galaxies; this is widely held to provide a better description of the halo statistics at high redshifts compared 
to the canonical Press & Schechter distribution. 

The result (see right-bottom panel in fig. 7) is to further reduce the excess to a factor ^2.5—3 for luminosities fainter 
than M1700 ~ —19. At present, two explanations may be offered to account for the residual excess. 

First, incompleteness in the data. This has been estimated by Vanzella et al. (2001) in the HDF-S data; at the faintest 
limits used to compute the LF (Iab — 27.25), it has been found to be close to 1.6 for extended sources. Even applying 
this correction factor to the whole faintest bin of the observed LF, the excess still remains. 

Second, additional sources of feedback; e.g., the feedback due to the photoionization of the intergalactic medium by 
photons escaping with high efficiency from stars and quasars (Benson et al. 2001) could reconcile the predictions with the 
observed LFs. In this context, the contribution of binary aggregations reduces the amount of feedback required to yield 
flat LFs as to match the observations. 

6. DISCUSSION, CONCLUSIONS, AND PERSPECTIVES 

We have developed a semi-analytic model (SAM) for galaxy evolution, which includes the effects of binary aggregations 
between satellite galaxies orbiting inside common DM haloes (see fig. 1). The corresponding evolution of the galaxy mass 
function (see the velocity functions in fig. 3) has been calculated in a A-CDM cosmology by solving in each host halo a 
kinetic Smoluchowski equation grafted onto our SAM code. 

The main effects of aggregations on the galaxy mass distribution are summarized as follows: 

i) Starting at z m 3 (see fig. 4), aggregations of satellite galaxies gradually deplete the number of low/intermediate mass 
galaxies, and flatten the slope of the mass function (see fig. 3). The aggregation cross section strongly depends on the 
ratio v/V re i between the galaxy circular velocity and the velocity dispersion of the halo wherein it orbits, so the increase 
of the average V re i following from the hierarchical clustering results in a progressive shift of the range of circular velocities 
where the aggregations flatten the mass distribution. At z = the resulting circular velocity distribution is flatter for 
small/intermediate circular velocities (50 <J v <; 200 km/s, see fig. 3). The result agrees with the outcomes from N-body 
simulations. 

ii) The total (dynamical friction plus binary aggregations) net merging rate peaks at z ~ 3.5 and declines sharply 
afterwards, in agreement with N-body results and consistent with observations. The contribution of satellite aggregations 
to the destruction rate (which mainly affects the low mass end of the mass distribution) becomes comparable to that 
produced by dynamical friction for z < 1, thus sustaining the galaxy destruction rate at low rcdhsifts. We have checked 
that the construction and destruction rates corresponding to the dynamics in our model (see fig. 4) agree with those 
obtained from N-body simulations within their limited mass resolution. 

The flattening of the mass distribution at small/intermediate masses velocities 50 < v <J 150 km/s) affects the galaxy 
luminosity functions as shown in figs. 6, 7. These, obtained on adopting for the stellar section the standard prescriptions 
of SAMs, have been compared with data from deep surveys with spectroscopic or photometric redshifts. The results 
concerning the luminosity functions are summarized as follows: 

• at low redshifts, the binary aggregations flatten the LF at small/intermediate luminosities, producing also an upturn 
at very faint magnitudes Mb > —16. At magnitudes brighter than L* the LFs are not appreciably affected. 

• at higher redshifts, the excess of the standard SAM predictions over the observed LFs at the faint end (by a factor 
<~ 8 at z ss 3 for Mi 700 ~ —18, sec Poli et al. 2001) is considerably reduced by the introduction of satellite aggregations 
computed in our present model. At z <~ 3, our LFs computed adopting the Sheth & Tormen (1999) mass distribution for 
the host DM haloes further reduces the excess by a factor w 3 in the faintest bins (—19.5 <J Mn 00 <; —18). 

The origin of the residual excess at high redshifts requires further investigation. Incompleteness at the faintest limits of 
the observed LF (Iab = 27.25) has been found to be around 1.6 for point and extended sources (Vanzella et al 2001). With 
regard to curing the residual excess, we note that changing the parameters in the canonical stellar feedback prescriptions 
is not a viable solution, since a larger feedback parameter (as required to flatten the galaxy LF still more) would result 
in a Tully-Fisher relation that is too faint when compared to the data. An alternative may be provided by sources of 
feedback other than SNe, e.g., the feedback due to photoionization of the intergalactic medium by stars and quasars (as 
advocated by Benson et al. 2001) could reconcile the predicted with the observed LFs. In this context, the contribution 
of binary aggregations to the flattening of the LF will alleviate the need for extreme feedback efficiencies. Further insight 
could be provided by the comparison of the model predictions with K-band luminosity functions and with the baryonic 
mass function of galaxies (as suggested by Salucci & Persic 1999), a point we will address elsewhere. 

On the other hand, our conclusions about the faint end of the galaxy LFs are little affected by other processes recently 
implemented in the SAMs, namely, the additional starbursts associated with the mergers following satellite aggregations 
(Somerville & Primack 1999), and the effect of tidal stripping of galactic subhaloes inside the host DM haloes (Benson et 
al. 2001). The former authors associate starbursts (with tunable timescale and amplitude) with each newly formed merger 
resulting from binary aggregations. We did not implement them in our model in order to single out dynamical effects of 
our description without introducing new (and uncertain) free parameters in the model. On the other hand, the starbursts 
associated with the mergers resulting from the galaxy encounters are expected to have a only a moderate inpact on the 
LF at faint /intermediate luminosities. In fact, the dynamics described by eqs. (3) and (4) is such that the aggregation 
events mainly occur between intermediate and small mass galaxies (the first being favoured by their larger cross section 
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and the latter by their large number), to form larger units; thus, the destruction term dominates at the small/intermediate 
masses while the construction term dominates at large masses. The starbursts associated with mergers are then expected 
to brighten mainly the bright end of the LF. 

As for tidal stripping, we note that Benson et al. (2001) have shown that it plays a minor role in determining the 
shape of the global galaxy luminosity function. This process, if anything, affects mainly the mass distribution at the very 
low-mass end (circular velocities v w 20 — 50 km/s), while satellite aggregations modelled here affect the mass distribution 
at larger scales (v s=s 50 — 150 km/s). 

We stress that our description of the dynamics of galaxies inside common DM haloes does not introduce new free 
parameters in the SAMs; so, it constitutes a step toward improving the physical treatement of the processes involved in 
galaxy formation and evolution, rather than merely enlarging the set of phenomenogical scaling laws and the associated list 
of free paremeters. A further step in such a direction will be constituted by a refined description of how central dominant 
galaxies form. In fact, although in the present paper the formation is treated in terms of dynamical friction to retain 
continuity with other works, in principle it could be included in the statistical description provided by the Smoluchowski 
kinetic equation. The process is due to the rare slow encounters for which the focussing term v/V re i contained in eq. (4) 
dominates. In this case, the V^-dependent cross section (i.e., eq. 4 without the velocity average) would become more 
than linear in mass in the low velocity tail of g(V re i) (in particular, X oc m 4//3 holds in such a case). In these conditions 
the aggregation dynamics described by eq. (3) are marked by a phase transition of gravitational nature which breaks 
the system of colliding galaxies, originally comprised of a number of comparable objects, into two phases: a prompt, 
dominant merger; and a number of satellites of much smaller mass, which gradually aggregate to the dominant object 
(see Cavaliere, Colafrancesco & Menci 1991; Menci, Colafrancesco & Biferale 1993). The transition takes place at time 
around Tdyn 

(m/M) 5 / 6 (the detailed expression is given by the above authors), which turns out to be comparable with 
the timescale actually adopted in the present paper, see eq. (2). 

Thus, in principle, the entire dynamics of galaxies in common haloes, including satellite aggregations and formation of 
central dominant mergers, could be treated within SAMs using only the kinetic description given in §3.2. In practice, how- 
ever, this would require the detailed consideration of the galaxy velocity distribution to describe properly the probability 
for the occurence of slow, resonant encounters leading to the phase transition; this goes beyond the average description 
provided by eq. (4). In the present paper the formation of a central merger is described only in terms of dynamical 
friction; this can be viewed as a mean field, average rendition of the critical phenomenon leading to the phase transition. 
In perspective, we plan to describe all galaxy interactions in terms of the kinetic theory; this will be the subject of a 
future paper. 

We gratefully acknowledge constructive advice of our referee toward correcting our manuscript and improving our 
presentation. Work supported by partial grants from ASI and MIUR. 
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APPENDIX A 



The statistics of the virialized DM condensations in the framework of the hierarchical clustering is commonly described 
in terms of the Extended Press & Schechter Theory (EPST; Bower 1991; Bond et al. 1991; Lacey & Cole 1993). Here we 
summarize the results used in the text. 

The number density of virialized structures of mass M at the cosmic time t is given by the Press & Shechter (1974) 
formula 



N H (M,t) = 



2 Po_ 
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where po is the cosmic average matter density, a(M) is the rms density fluctuations of the linear perturbation field in 
spheres of mass M (see, e.g., Lacey & Cole 1993), and S c (t) is the density threshold for collapse. For a(M) we adopt the 
form corresponding to Cold Dark Matter (CDM) density perturbations in a flat Universe dominated by a cosmological 
constant Cl\ = 0.7. The threshold 5 c (t) corresponds to the value - extrapolated to the present using the linear growth 
factor I?(t, Oo,^a) for the density perturbations - of the overdensity of a homogeneous sphere at the point where the 
exact non-linear theory predicts collapse to a singularity. Its normalization 5 c q at z = and its time evolution I? _1 (i) 
depend on cosmology; in a critical = 1 Universe S c0 = 1.686 and D^ 1 (z) = (1 + z) hold; the corresponding expressions 
for a flat Universe with a nonzero cosmological constant are given by, e.g., Eke, Cole & Frenk (1996). 

Since the variance of the density field is an inverse function of M and the density threshold S c (t) lowers with time, 
larger and larger overdensities collapse to form structures of increasing mass (hierarchical clustering). For a given DM 
mass M it is possible to compute the merging rates and the progenitor distributions. 

In particular, the probability that a condensation of mass M present at time t form a halo of mass between M' and 
M' + dM' at time t + dt is given by the following expression: 
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where a and a' are the rms density fluctuations corresponding to the masse M and M' , respectively. 

The history of the previous mergers experienced by a mass M is instead described by the probability distribution that 
a given mass M at time t has a progenitor of mass M' at time t' < t: 
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The survival time of a halo of mass M is defined as the cosmic time at which its mass has grown to qM, and hence 
depends on the choice of the "mass step" of the merging tree. For a generic step parameter q, the probability prob[ri < r] 
that a halo with mass M at time t has a survival time t + n with t; < r reads (Lacey & Cole 1993): 



probWi < 



1 Mt) 
2 



S c (t) 
x [1 - erf(X)] - 



2SJt + T)] !M±I] feW-^c(t+r)] 



[1 - erf(Y)] 



o 2 {gM) [6 c (t) - 2 S c (t + r)] + a 2 (M) S c (t + r) 
{2a 2 (M) a 2 (qM) [a 2 (M)-a 2 (qM)}} 1 / 2 
a 2 (M)d c (t + r)-a 2 (qM)S c (t) 
{2 a 2 {M) a 2 {qM) [a 2 {M) - a 2 (qM)}y/ 2 



(A-4a) 
(A-4b) 

(A-4c) 



To make contact with previous works (Cole et al. 1994; Cole et al. 2000), our default choice for the mass step is q = 2. 
We have verified that the SAM results are robust to smaller values down to q = 1.2, as discussed by Cole et al. (2000). 

Finally, we recall the expression for the mass function recently proposed by Sheth & Tormen (1999; see also Sheth, Mo 
& Tormen 2001) to provide a better fit to the N-body results. This has the form 



N H (M,t) = A 




dlna 



dlnM 



aS 2 (t) 



SJt) - "*| (t) 



a{M) 



(A-5) 



The values of the fitting coefficients A = 0.3222, a — 0.707 and p = 0.3 are obtained by the above authors from comparison 
with the N-body results. The agreement of the above expression for Nh{M) with the simulations has been confirmed by 
an independent analysis (Jenkins et al. 2001). 
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APPENDIX B 



We derive the average value for the tidal radius rud of galaxies orbiting inside a host halo with radius R and circular 
velocity V. Such a radius is that appropriate for a galactic subhalo which survives the tidal stripping of the host halo, a 
condition which requires the density of the galactic subhalo within rud to exceed the density of the host halo interior to 
the pericentre r p of its orbit. 

To this end we adopt the approach of Ghigna et al. (1998) who showed how the condition for survival against tidal 
stripping translates approximately into r t id ~ r p{ v /V)- The above expression for r t %d has been tested by the above 
authors against high resolution N-body simulations, and has been proved to agree with the values of rtid measured in the 
simulations for all subhaloes except the few on very eccentric orbits (the latter have measured radii larger than expected 
partly due to the formation of tidal tails). 

We average the above relation for r t id over the distribution of pericentres obtained by Ghigna et al. (1998) from N-body 
simulations, which we fit with a modified lognormal expression (in the variable (r p /R — 0.08) with logarithmic mean -1.3 
and logarithmic variance 0.6). Performing the average yields for r U d the following expression: 



where the lower limit r cut corresponds to the minimum tidal radius that a galactic subhalo can have without being severely 
distorted or disrupted. Following Bullock, Kravstov and Weinberg (2000), we adopt for r cut the radius of the peak of 
the galaxy velocity profile. For a Navarro et al. (1997) circular velocity profile r cut — 2.16r 2 oo/c holds, where c is the 
concentration parameter of the subhalo, and f2oo is the radius where the average density of the subhalo would equal 
200 p c . The equations relating c and r 2 oo with the circular velocity v of the subhalo are given by the above authors. The 
validity of the above value for r cut has been tested against N-body results by Bullock et al. (2000). 

It is easy to recast the relation rud — aRv/V in terms of density. Substituting the relation v/V = (m/M) 1 / 2 (R/rud) 1 ^ 2 
one obtains (r t id/R) 3 = a 2 m/M so that the ratio of the subhalo to the host halo average densities is given by I /a 2 . The 
average over the distribution of pericentres defining the value of a yields typical values for such a ratio ranging between 
3 and 5, depending on the concentration parameter c of the subhaloes. 

Note that taking the tidal radius r t id — r p v/V as the limiting radius of the subhaloes is fully appropriate only for 
subhaloes much smaller than their host halo. This is the case for the bulk of the population of satellite galaxies (the one 
mainly contributing to binary aggregations) for which our treatment constitutes a valid approximation. 




(B-l) 



Menci et al. 
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